# en caso no se tengan las dependencias
# install.packages("devtools")
# devtools::install_github("r4ss/r4ss", ref="development")
# install.packages("caTools")
# library("caTools")
# # install.packages("r4ss")
library(r4ss)
library(here)
#remotes::install_github("PIFSCstockassessments/ss3diags")
library(ss3diags) # diagnosticos del modelo
library(ss3sim) # evaluación de sesgo
library(kableExtra)
library(tidyverse)Implementación de evaluación de stock en reineta Brama australis en SS3
ANTECEDENTES
Descripción y objetivo del documento
El principal objetivo de este documento es presentar los avances en la implementación del modelo de evaluación integrado para reineta Brama australis en la zona centro sur de Chile. Este proceso de evaluación ha sido implententado en Stock Synthesis con la finalidad de integrar distintas piezas de información disponible, analiar desempeño de los modelos y escenarios probados y por último, avanzar en la discusión sobre la utilidad de este tipo de aproximaciones metodologicas y vinculación con la asesoria de IFOP hacia el manejo.
Por lo tanto, este documento contiene el flujo de análisis y modelación de los distintos escenarios de Brama australis Reineta como parte de la asesoría técnica que lleva a cabo el IFOP, mediante el Grupo Técnico Demersal y Crustáceos del Departamento de Evaluación de Recursos.
METODOLOGÍA
Área de estudio
El área de estudio comprende la principal área de operación de la flota arrastrera correspondiente a la zona centro sur de Chile entre el límite norte de la IV región y los 41°28’S. Donde a lo largo de la costa operan las flotas artesanal con redes de enmalle y espinel, mientras que la flota de industrial desarrolla su operación de captura con arrastre.
Modelo conceptual
El modelo de la dinámica poblacional de reineta, corresponde a un enfoque de evaluación del tipo estadístico con estructura de edad, donde la dinámica progresa avanzando en el tiempo t, y las capturas son causantes de la mortalidad por pesca F, la mortalidad natural es constante M = 0,35. La relación entre la población y las capturas responde a la base de la ecuación de Baranov, y se consideran para el modelo y estimaciones el rango de edad entre 2 a 12+ (años). Sin embargo, las estimaciones del modelo tienen su origen en la edad cero sobre la base de una condición inicial estado estable. La dinámica esta modelada por un reclutamiento tipo Beverton y Holt.
Plataforma de modelación
Los modelos implementados fueron configurados utilizando Stock Synthesis (SS3 de aquí en mas)(https://vlab.noaa.gov/web/stock-synthesis) con la versión mas actualizada (V.3.30.18). SS3 es un modelo de evaluación de stock edad y talla estrucuturado, en la clase de modelos denominados “Modelo de análisis integrado de evaluación de stock”. SS3 tiene un sub-modelo poblacional de stock que simula crecimiento, madurez, fecundidad, reclutamiento, movimiento, y procesos de mortalidad, y sub-modelos de observation y valores esperados para diferentes tipos de datos. El modelo es codificado en C++ con parámetros de estimación activados por diferenciación automática (ADMB) (Fournier et al., 2012; Methot & Wetzel, 2013). El análisis de resultados y salidas emplea herramientas de R e interfase gráfica de la librería r4ss (https://github.com/r4ss/r4ss) (Taylor, 2019).
Parámetros biológicos
La Tabla identifica los parametros de información disponbles que provienen del reporte técnico de Contreras (2020). En este informe no se realizó una revisión de los parámetros biologicos disponibles, dado que se utilizaron los de las evaluaciones previas para ser consistentes con los informes anteeriores de evaluación y la información utilizada.
| LO | HI | INIT | PRIOR | PR_SD | PR_type | PHASE | env_var&link | |
|---|---|---|---|---|---|---|---|---|
| NatM_p_1_Fem_GP_1 | 5e-02 | 0.900000 | 0.350 | 0.350 | 0.60 | 3 | -4 | 0 |
| L_at_Amin_Fem_GP_1 | 1e+00 | 35.000000 | 13.000 | 13.000 | 0.50 | 0 | -1 | 0 |
| L_at_Amax_Fem_GP_1 | 2e+01 | 80.000000 | 49.000 | 49.000 | 99.00 | 0 | -1 | 0 |
| VonBert_K_Fem_GP_1 | 8e-02 | 0.550000 | 0.180 | 0.180 | 99.00 | 0 | -1 | 0 |
| CV_young_Fem_GP_1 | 1e-01 | 0.600000 | 0.300 | 0.300 | 99.00 | 0 | -2 | 0 |
| CV_old_Fem_GP_1 | 1e-01 | 0.400000 | 0.150 | 0.200 | 99.00 | 0 | -2 | 0 |
| Wtlen_1_Fem_GP_1 | -3e+00 | 3.000000 | 0.021 | 0.021 | 99.00 | 0 | -5 | 0 |
| Wtlen_2_Fem_GP_1 | -3e+00 | 6.000000 | 2.940 | 2.940 | 99.00 | 0 | -5 | 0 |
| Mat50%_Fem_GP_1 | 2e+01 | 54.000000 | 37.000 | 37.000 | 99.00 | 0 | -5 | 0 |
| Mat_slope_Fem_GP_1 | -3e+00 | 3.000000 | -0.200 | -0.200 | 99.00 | 0 | -5 | 0 |
| Eggs_intercept_Fem_GP_1 | -3e+00 | 3.000000 | 1.000 | 1.000 | 0.05 | 0 | -50 | 0 |
| Eggs_slope_wt_Fem_GP_1 | -3e+00 | 3.000000 | 0.000 | 0.000 | 0.05 | 0 | -50 | 0 |
| CohortGrowDev | 1e-01 | 10.000000 | 1.000 | 1.000 | 0.05 | 0 | -1 | 0 |
| FracFemale_GP_1 | 1e-06 | 0.999999 | 0.500 | 0.500 | 0.05 | 0 | -99 | 0 |
El reclutamiento fue modelado mediante una curva logística de Beverton y Holt, con los parámetros indicados en la Tabla .
step <- ctl$SR_parms[1:2,1:7]
kbl(step, booktabs = T,format = "html",
caption = "\\label{t2}Parámetros Relación S-R") %>%
kable_styling(latex_options = "HOLD_position")| LO | HI | INIT | PRIOR | PR_SD | PR_type | PHASE | |
|---|---|---|---|---|---|---|---|
| SR_LN(R0) | 3.0 | 31 | 8.815050 | 10.3 | 10.00 | 0 | 1 |
| SR_BH_steep | 0.2 | 1 | 0.614248 | 0.7 | 0.05 | 1 | 4 |
Datos utilizados
Desembarque industrial y artesanal del período (1994-2021) separados por flota, provenientes de las estadisticas oficiales de Sernapesca (Subsecretaria de Pesca, 2021). Al disponer de los desembarques oficiales por flota, es posible segregar información oficial por flotas, siendo factible a la vez disponer de datos oficiales (reportados). Cabe señalar que en esta pesquería aun no se realizan procesos de corrección de desembarques.
Información del Programa de Seguimiento de la pesquería de pesquerías demersales del Instituto de Fomento Pesquero.
La información proviene del monitoreo artesanal e industrial en la zona centro-sur de Chile, en donde se destacan dos flotas de pesca, la artesanal de enmalle y artesanal de espinel, siendo esta ultima la mas importante en terminos de registros e historial.
En ambos casos se obtienen datos de: i) estructura de tamaños, ii) composiciones por edad, iii) parámetros de crecimiento y iv) peso anuales por edad/talla y años.
Por otro lado, los rendimientos de pesca de cada flota fueron estandarizados mediante modelos lineales generalizados.
- En el caso de la pesquería industrial, el monitoreo permitió obtener composición de abundancia a la edad entre los años 2017 y 2021.
El esquema de evaluación presentado considera una modelación secuencial por flotas artesanal e industrial. La flota con mayor historial pesquero de acuerdo a los analisis es la flota artesanal de espinel, que para este caso se determinara como modelo base de evaluuación.
La evaluación de stock en SS3 de Brama australis en la zona Centro Sur de Chile se realiza de manera jerarquica integrando la información relativa a las tres flotas que operan en el recurso. A saber;
- Espinel Artesanal (1)
- Enmalle Artesanal (2)
- Industrial (3)
En la Tabla 3. se enumeran y describen los escenarios de modelación en función de las flotas disponibles.
| Escenario | Descripción |
|---|---|
| s1 | Todas las flotas |
| s1a | Todas las flotas (comps Edades Industrial) |
| s2 | Flota Artesanal Enmalle |
| s3 | Flotas artesanales (Espinel y Enmalle) |
| s4 | Flota Artesanal Espinel (Modelo Base) |
| s5 | Indice Biomasa Zhou (Zhou et al., 2009) |
Para avanzar en la implenteación metodológica, se establece con fines comparativos modelo por flotas artesanales, donde un modelo utiliza la información de enmalle artesanal, para luego sumar la flota enmalle artesanal, para terminar incorporando la información de la flota industrial. Tambmien se ehecutan las flotas por separado (Tabla 1).
Modelos reportados
Cabe mencionar que el modelo utilizado para la toma de decisiones esta basado en un modelo de producción basado en datos de captura con uso de supuestos asociados a los niveles de agotamiento (Contreras, 2020; Zhou et al., 2009). Por otro lado, se ha avanzado en la implenteación de modelos integrados utilizando datos de la flota artesanal de enmalle (Contreras, 2020).
Al momento de la elaboración de este informe, existen tres modelos implementados a la fecha, a saber el s2, s3 y s4. Estos modelos seran presentados a continuación en referencia las piezas de información disponibles y salidas de la evaluación.
Datos utilizados en s2
El modelo s2 contiene los datos de enmalle de la flota artesanal de reineta que se presentan a continuación en la Figura 2.
SSplotData(base.model2, subplot = 1,
fleetnames = c("Enmalle"),
fleetcol = c(5))Desembarques asociados a las flotas del modelo s2.
SSplotCatch(base.model2, subplots = 2,
fleetnames = c( "Enmalle"),
fleetcols = c(5), forecastplot = T)Datos utilizados en s3
El modelo s2 contiene los datos de enmalle y espinel de la flota artesanal de reineta que se presentan a continuación en la Figura 2.
Desembarques asociados a las flotas del modelo s3.
SSplotCatch(base.model3, subplots = 2,
fleetnames = c("Espinel", "Enmalle"),
fleetcols = c(2,5), forecastplot = T)Datos utilizados en s4
El modelo s4 contiene los datos de espinel de la flota artesanal de reineta que se presentan a continuación en la Figura 4.
SSplotData(base.model4, subplot = 1,
fleetnames = "Espinel",
fleetcol = 2)Desembarques asociados a las flotas del modelo s4.
SSplotCatch(base.model4, subplots = 2,
fleetnames = c("Espinel"),
fleetcols = 2)CPUE
Las tablas 3, 4 y 5 muestran las series de rendimientos nominales definidas por las pesquerías de Enmalle y Espinel para reineta.
| year | seas | index | obs | se_log |
|---|---|---|---|---|
| 1998 | 6 | 1 | 667 | 0.3 |
| 1999 | 6 | 1 | 1019 | 0.3 |
| 2000 | 6 | 1 | 1108 | 0.3 |
| 2001 | 6 | 1 | 1142 | 0.3 |
| 2002 | 6 | 1 | 530 | 0.3 |
| 2003 | 6 | 1 | 297 | 0.3 |
| 2004 | 6 | 1 | 360 | 0.3 |
| 2005 | 6 | 1 | 732 | 0.3 |
| 2006 | 6 | 1 | 590 | 0.3 |
| 2007 | 6 | 1 | 725 | 0.3 |
| 2008 | 6 | 1 | 654 | 0.3 |
| 2009 | 6 | 1 | 427 | 0.3 |
| 2010 | 6 | 1 | 1252 | 0.3 |
| 2011 | 6 | 1 | 1847 | 0.3 |
| 2012 | 6 | 1 | 1437 | 0.3 |
| 2013 | 6 | 1 | 672 | 0.3 |
| 2014 | 6 | 1 | 1006 | 0.3 |
| 2015 | 6 | 1 | 694 | 0.3 |
| 2016 | 6 | 1 | 402 | 0.3 |
| 2017 | 6 | 1 | 627 | 0.3 |
| 2018 | 6 | 1 | 654 | 0.3 |
| 2019 | 6 | 1 | 1027 | 0.3 |
| 2020 | 6 | 1 | 1265 | 0.3 |
| 2021 | 6 | 1 | 1265 | 0.3 |
| year | seas | index | obs | se_log | |
|---|---|---|---|---|---|
| 25 | 2006 | 6 | 2 | 36.26 | 0.3 |
| 26 | 2007 | 6 | 2 | 97.54 | 0.3 |
| 27 | 2008 | 6 | 2 | 565.45 | 0.3 |
| 28 | 2009 | 6 | 2 | 306.94 | 0.3 |
| 29 | 2010 | 6 | 2 | 155.50 | 0.3 |
| 30 | 2011 | 6 | 2 | 858.38 | 0.3 |
| 31 | 2012 | 6 | 2 | 513.68 | 0.3 |
| 32 | 2013 | 6 | 2 | 260.25 | 0.3 |
| 33 | 2014 | 6 | 2 | 353.84 | 0.3 |
| 34 | 2015 | 6 | 2 | 226.97 | 0.3 |
| 35 | 2016 | 6 | 2 | 181.68 | 0.3 |
| 36 | 2017 | 6 | 2 | 132.54 | 0.3 |
| 37 | 2018 | 6 | 2 | 169.04 | 0.3 |
| 38 | 2019 | 6 | 2 | 276.58 | 0.3 |
| 39 | 2020 | 6 | 2 | 226.65 | 0.3 |
| 40 | 2021 | 6 | 2 | 741.35 | 0.3 |
| year | seas | index | obs | se_log | |
|---|---|---|---|---|---|
| 41 | 2011 | 7 | 3 | 584 | 0.3 |
| 42 | 2012 | 7 | 3 | 430 | 0.3 |
| 43 | 2013 | 7 | 3 | 139 | 0.3 |
| 44 | 2014 | 7 | 3 | 713 | 0.3 |
| 45 | 2015 | 7 | 3 | 528 | 0.3 |
| 46 | 2016 | 7 | 3 | 320 | 0.3 |
| 47 | 2017 | 7 | 3 | 319 | 0.3 |
| 48 | 2018 | 7 | 3 | 395 | 0.3 |
| 49 | 2019 | 7 | 3 | 643 | 0.3 |
| 50 | 2020 | 7 | 3 | 707 | 0.3 |
| 51 | 2021 | 7 | 3 | 729 | 0.3 |
Por otro lado, se presentan los graficos que señalan las tendiencias de los rendimientos (Figura 5).
RESULTADOS
Respecto a los valores y parámetros biológicos modelados, los siguientes gráficos identifican los estimadores puntuales del recurso en s2, s3 y s4. En el primer analisis visualizamos el modelo de crecimkento individual del recurso en cada modelo.
SSplotBiology(base.model2, subplots =1, labels = c("Length (cm)", "Age (yr)", "Maturity", "Mean weight (kg) in last year",
"Spawning output", "Length (cm, beginning of the year)", "Natural mortality",
"Female weight (kg)", "Female length (cm)", "Fecundity", "Default fecundity label",
"Year", "Hermaphroditism transition rate", "Fraction females by age at equilibrium"),
)SSplotBiology(base.model3, subplots =1, labels = c("Length (cm)", "Age (yr)", "Maturity", "Mean weight (kg) in last year",
"Spawning output", "Length (cm, beginning of the year)", "Natural mortality",
"Female weight (kg)", "Female length (cm)", "Fecundity", "Default fecundity label",
"Year", "Hermaphroditism transition rate", "Fraction females by age at equilibrium"),
)SSplotBiology(base.model4, subplots =1, labels = c("Length (cm)", "Age (yr)", "Maturity", "Mean weight (kg) in last year",
"Spawning output", "Length (cm, beginning of the year)", "Natural mortality",
"Female weight (kg)", "Female length (cm)", "Fecundity", "Default fecundity label",
"Year", "Hermaphroditism transition rate", "Fraction females by age at equilibrium"),
)Ajustes
Una de las mas básicas formas de identificar el desempeño de los modelos es mirando los ajustes asociados a la estimación de parámetros y variables. En primer lugar observamos los ajustes a las tallas de cada escenario.
SSplotComps(base.model2, subplots = 1,
fleetnames = c("Enmalle"),
sizemethod = 2,
smooth = TRUE)SSplotComps(base.model3, subplots = 1,
fleetnames = c("Espinel", "Enmalle"),
sizemethod = 2,
smooth = TRUE)SSplotComps(base.model4, subplots = 1,
fleetnames = c("Espinel"),
sizemethod = 2,
smooth = TRUE)Biomasa desovante
La estimación de biomasa total por modelo da cuenta de las diferencias entre modelos. A continuación se grafican las salida de las biomasas de los modelos s2, s3 y s4.
SSplotTimeseries(base.model2, subplot = 1)SSplotTimeseries(base.model3, subplot = 1)SSplotTimeseries(base.model4, subplot = 1)Reclutamiento
Las siguientes figuyras muestran los desvios del reclutamiento para los escenarios s2, s3 y s4.
SSplotRecdevs(base.model2, subplot = 2)SSplotRecdevs(base.model3, subplot = 2)SSplotRecdevs(base.model4, subplot = 2)Mortalidad por pesca
La mortalidad por pesca se identifica en las siguientes figuras.
SSplotSummaryF(base.model2)SSplotSummaryF(base.model3)SSplotSummaryF(base.model4)Retrospectivo
Los análisis retrospectivo, dan cuenta de diferencias en los patrones entre modelo 1 y 2. En el caso del reclutamiento el patrón corresponde a una sub-estimación (Figura 14) en M1, mientras que en M2 el patrón se relaciona con una sobre-estimación (Figura 15).
mydir <- dir2
SS_doRetro(
masterdir = mydir,
oldsubdir = "",
newsubdir = "retrospectives",
years = 0:-5
)

mydir3 <- dir3
SS_doRetro(
masterdir = mydir3,
oldsubdir = "",
newsubdir = "retrospectives",
years = 0:-5
)

mydir4 <- dir4
SS_doRetro(
masterdir = mydir4,
oldsubdir = "",
newsubdir = "retrospectives",
years = 0:-5
)

Comparación
La siguiente Tabla muestra los componentes de la probabilidad asociados a la estimación de cada escenario testeado en este documento. A su vez se identifican los parámetros estimados por cada assessment. Tambien podemos identificar las diferencias entre modelos en las principales variables poblacionales estimadas para lostres escenarios (s2, s3 y *s4)
mod.sum <- SSsummarize(list(base.model2, base.model3, base.model4))SSplotComparisons(mod.sum,
legendlabels=c("s2", "s3", "s4"),
print = TRUE,
plot=TRUE,
png = TRUE, models = "all",
plotdir = diri)kbl(comp, booktabs = T,format = "html",
caption = "Comparacion likelihood y parámetros s2, s3 y s4") %>%
kable_styling(latex_options = "HOLD_position")| Label | s2 | s3 | s4 |
|---|---|---|---|
| TOTAL_like | 273.363000 | 3.99920e+03 | 774.834000 |
| Survey_like | 26.261800 | 2.28558e+03 | 8.293300 |
| Length_comp_like | 244.007000 | 1.66479e+03 | 713.349000 |
| Parm_priors_like | 2.862110 | 2.78771e+00 | 2.992780 |
| Recr_Virgin_thousands | 59.431200 | 1.86219e+08 | 20.085500 |
| SR_LN(R0) | 4.084820 | 1.90424e+01 | 3.000000 |
| SR_BH_steep | 0.987895 | 9.42771e-01 | 0.999218 |
| NatM_uniform_Fem_GP_1 | 0.350000 | 3.50000e-01 | 0.350000 |
| L_at_Amax_Fem_GP_1 | 49.000000 | 4.90000e+01 | 49.000000 |
| VonBert_K_Fem_GP_1 | 0.180000 | 1.80000e-01 | 0.180000 |
| SSB_Virgin_thousand_mt | 22.023000 | 6.90054e+07 | 7.443000 |
| Bratio_2021 | 2.369070 | 1.92201e+00 | 0.224274 |
| SPRratio_2021 | 0.636230 | 1.20000e-06 | 0.857473 |
DISCUSIÓN
Asociadas a la implementación metodológica
La implementación metodologica de la evaluación de stock mediante un modelo integrado de reineta para la zona Centro Sur de chile presenta multiples desafíos, y que algunos componentes de este proyecto fueron atendidos en este documento. El uso de plataformas de evaluación de stock como SS3 (Methot & Wetzel, 2013) ha permitido explorar otros escanarios de la naturaleza asociada a este recurso, dinámica poblacional y pesquería, en este caso, incorporar nuevas fuentes y piezas de datos disponibles los cuales dieron paso a los resultados presentados.
Datos y piezas de información disponibles
Este trabajo de implementación metodológica fue posible dado el trabajo de análisis de datos disponibles de las tres flotas que operan sobre este recurso, a saber: espinel, enmalle e industrial. Este trabajo ha sido desarrollado por el Grupo de Trabajo de Demersales del Departamento de Evaluación de Recursos de IFOP. Los principales avances tienen relación con obtención de rendimientos estandarizados y estructuras de tallas asociadas a cada flota. Estos datos presentan un avance en la información disponible para la implentación de un modelo integrado con dinamica en edades como se plantea en este documento.
Asociadas a la evaluación de stock
Respecto a los resultados de la evaluación de stock, es posible identificar la consistencia en la estimación de los modelos s2 y s3, es decir, los modelos que utilizaron las flotas de enmalle y espinel por separado respectivamente. Por otro lado, el modelo que combinó las flotas de enmalle y espinel (s3) fue el que tuvo el desempeño mas bajo de los tres modelos presentados.
Trabajo en progreso
Este trabajo presenta componentes que aun deben ser discutidos y analizados de manera mas extensa. En primer lugar, es necesario parametrizar en el ambito biologico de la especie y de la pesquería los tres modelos analizados. En segundo lugar, se deben terminar de implementar los modelos con la información de la flota industrial, los cuales deben ser integrados a la información disponible de las otras flotas artesanales (modelo s1 y s1a). Por último, se debe avanzar en un modelo que utilice otro indicador de abundancia alternativo al de la CPUE. En este caso hemos propuesto utilizar las salidas del modelo de producción utilizado para el manejo y que actualmente se presenta en la evaluación (Zhou et al., 2009)
Consideramos que este trabajo se sitúa en el camino correcto para la integración de las piezas de información disponible en esta pesquería, sin embargo, la implenteación y mejoras presentadas en este documento deben ser analizadas a la luz de un trabajo que no ha terminado y en progreso.